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Abstract. Considering a gas of self-propelled particles with binary interactions, 
we derive the hydrodynamic equations governing the density and velocity fields 
from the microscopic dynamics, in the framework of the associated Boltzmann 
equation. Explicit expressions for the transport coefficients are given, as a function 
of the microscopic parameters of the model. We show that the homogeneous 
state with zero hydrodynamic velocity is unstable above a critical density (which 
depends on the microscopic parameters), signaling the onset of a collective motion. 
Comparison with numerical simulations on a standard model of self-propelled 
particles shows that the phase diagram we obtain is robust, in the sense that 
it depends only slightly on the precise definition of the model. While the 
homogeneous flow is found to be stable far from the transition line, it becomes 
unstable with respect to finite-wavelength perturbations close to the transition, 
implying a non trivial spatio-temporal structure for the resulting flow. We find 
solitary wave solutions of the hydrodynamic equations, quite similar to the stripes 
reported in direct numerical simulations of self-propelled particles. 



PACS numbers: 05.70.Ln, 05.20.Dd, 64.60.Cn 

1. Introduction 

In the recent years, a lot of effort has been expended with the aim of explaining the 
collective behaviour of living systems ]T] . Such collective behaviours can be observed 
on many different scales including mammal herds [2], crowds of pedestrians (3] H], 
bird flocks [5], fish schools [6], insect swarms [7], colonies of bacteria [8], molecular 
motors [10] and even interacting robots [IT] , It turns out that the collective 
properties of such systems seem to be quite robust and universal. Accordingly, this 
field attracted the interest of the statistical physics community with the challenge of 
introducing minimal models that could capture the emergence of collective behaviour. 
One important class of models consists of the so-called self-propelled particles models, 
for which the onset of collective motion without a leader is present. Vicsek et al. [T^inB] 
introduced a simple model defined on a continuous plane, where agents (or animals) 
are represented as point particles with a velocity of constant amplitude. Noisy 
interaction rules tend to align the velocity of any given particle with its neighbors. 
Extensive numerical simulations of this model have been performed |14[ 115] , showing 
the presence of a phase transition from a disordered state, at high enough noise, 
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to a state with collective motion. A different approach is to consider the problem 
at a coarsed-grained level and to describe the dynamics in terms of hydrodynamic 
fields. The equations governing the evolution of these hydrodynamic fields can be 
either postulated phenomenologically [16] , on the basis of symmetry and conservation 
laws considerations (TTl [18] , or derived from specific microscopic models [19l [20] . The 
equations of motion of the hydrodynamic field are derived from the microscopic model 
through a Boltzman approach. 

Following an earlier publication [19], the motivation of the present work is 
to derive, from a microscopic model, the hydrodynamic equations describing at a 
coarse-grained level the flow of self-propelled particles (SPP), and to compare the 
resulting description with numerical simulations of an agent-based model of SPP. The 
analytical framework we use is that of the Boltzmann equation. Accordingly, a suitable 
microscopic model for such a treatment is a continuous time model with interactions 
reducing to binary collisions. In order to show that the most salient features of the 
coarse-grained analytical description are not specific to a binary collision model, we 
use for the numerical simulations a standard agent-based model [HI [13], that has 
been well characterized in the litterature [TH [15] . Note that some comparison with 
numerical simulations of an agent-based model with binary interaction have already 
been presented in Ref. [19]. 

2. Microscopic models of interacting self-propelled particles 

2.1. Definition of the models 

2.1.1. Continuous time model with binary collisions. Following Ref. 19], we 
introduce a simple model that captures the essential physics of assemblies of self- 
propelled particles, while being suitable for a description in terms of a Boltzmann 
equation. We consider the evolution of self-propelled point-like particles on a two- 
dimensional plane. The displacement of each particle i is governed by a velocity vector 
Wi . In order to account for the self-propelling property, we assume that the modulus of 
the velocity vector is fixed to a value vo, identical for all the particles, so that only the 
direction of the vector plays a role in the dynamics. The relevant dynamical variables 
are then the angles 9i that the vectors form with a fixed reference direction. It 
is important to note at this stage, that fixing the modulus of the velocity breaks the 
Galilean invariance of the system. Hence one should not expect that the eventually 
obtained hydrodynamic equations obey such an invariance, contrary to what happens 
in usual flows. 

Apart from the ballistic evolution according to their velocity vector, particles also 
experience stochastic events that punctuate their dynamics. These stochastic events 
are of two different types. The simplest ones are self-diffusion events, that is, the angle 
9 of an isolated particle changes, with a probability A per unit time, to 9' = 9+r], where 
r\ is a noise with distribution po(rj) and variance Cq. In the following, we consider a 
Gaussian distribution for po(v)> a l so taking into account the periodicity of 9. This 
type of stochastic events lead to a diffusive behaviour at large scale, thus preventing 
the system from having a trivial (pseudo)collective motion, of purely ballistic nature. 
To drive the system into an organised state where a genuine collective motion sets in, 
one has to introduce interactions between the particles. Given that we wish to use a 
Boltzmann approach to study the model, it is natural to consider binary interactions 
between particles. These binary interactions are introduced as follows. When two 
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particles get closer than a threshold distance do, their velocity angle Q\ and 6* 2 are 
changed into 8[ and 9' 2 according to: 

0[ = 6 + ^21:1 e' 2 = 6 + m [2nl (1) 

where 8 is defined by: 

6 = arg(e i91 + e 102 ) . (2) 

The noises r\\ and r\i are independent Gaussian variables with variances a 2 . Note 
that cr 2 may differ from the variance a 2 of the noise associated to the self-diffusion of 
particles. 



2.1.2. Agent-based model for numerical simulations. In order to compare the results 
of the analytical approach based on the binary collision model to direct numerical 
simulations, we use a slight generalization of the standard Vicsek model [121 EES]- 
The motivation for simulating numerically a model different from the one we used in 
the analytical approach is twofold. First, the Vicsek model has been thoroughly 
characterized in the literature [131 HH EE], making it a useful benchmark for 
comparison. Second, and most importantly, a model with continuous time dynamics 
and binary collisions is well-suited for a Boltzmann equation approach, but very 
inefficient from the point of view of direct numerical simulations. In constrast, the 
Vicsek model, with a discrete time dynamics and multi-neighbour interactions, is much 
more efficient to simulate. 

The agent-based model we consider consists in TV particles on a two-dimensional 
space of area L x L, with periodic boundary conditions. Each particle j at any instant 
t has a constant modulus speed vq. This property allows the mapping of velocity on 
complex numbers. Then a particle is located by a two-coordinate vector x^- and an 
angle which gives its speed direction. We define the vicinity V* of j at time t as the 
disk centred on j with a radius do- Then the direction of j at the next instant t + At 
is simply the direction of the averaged speed over all particles which are embedded in 
its vicinity, including j itself, up to a noise term. If there is no neighbour in the disk 
of interaction, self-diffusion occurs randomly: 



nt + At 

j 



arg 



Efeevj^j +t$, ifVj^{j}, 
7 ? o£*, with probability AAt, if Vj = {j}, ( 3 ) 



qt with probability 1 - AAt, if V* = {j}, 



X i+At _ x t , „,_„^t+At^ 



, 3> " " " J 

v o e(0t +At )At, (4) 

where e(0) is the unit vector of direction 9. The parameters rj and 770 are the noise 
amplitudes for collision and self-diffusion respectively. The random number £j is 
uncorrelated in time and from one particle to another. Its distribution is flat on 
[— 7r,7r]. The slight generalization with respect to the standard Vicsek model consists 
in the introduction of the parameter A, which characterizes the probability of self- 
diffusion per unit time. In the original model, AAt = 1. Note that, whenever possible, 
we have defined the agent-based model with notations consistent with that of the 
binary collision model, in order to facilitate comparison between the two models. 

The Vicsek model has been studied in details in the literature [TH [T4J [15] . 
A transition toward collective motion has been reported in early studies |12j . and 
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later shown to exhibit strong finite size effects jT3] . In |Appendix A[ we recall the 
methodology used to study the transition, and in particular the finite size scaling 
effects. 

2.2. Dimensionless parameters 

A first step in the understanding of the models is to identify the relevant dimensionless 
parameters and the possible regimes. Let us first consider the different length scales 
appearing in this problem: the interaction range do, the ballistic length £bai = vq/X, 
and the typical distance between particles £ pp — l/^fp. With these three different 
lengths, one can form the following dimensionless numbers H and B: 

h= ^pp = 1 B= t _^L = ^L . ( 5 ) 

do d 0y /p ' d d X ' 

H characterises whether a system is diluted (H ^> I) or dense. One can see B as 
the relative weight of stand-alone flight over interaction. If B is large, ballistic flight 
is more important than collision and we can expect that particles are less correlated 
locally. 

These numbers turn out to play an important role in the identification of the 
regimes of validity of the approximations we use, as seen in the following. The model 
also exhibits different behaviours for the different regimes which are defined by these 
numbers. At fixed noise intensity and fixed B, a more (resp. less) dense system is 
expected to move (resp. not) in a collective manner. At fixed noise and fixed dilution 
H, increasing B makes the flights more ballistic, which should favor collective motion. 
So one can guess that a relevant control parameter will be a combination of H and B 
(see section |47T|) . 

2. 3. Summary of the main results 

The paper is organised as follows. Section [3] is devoted to the derivation from the 
Boltzmann equation, through a specific approximation scheme, of the hydrodynamic 
equations for the continuous time binary collision model. Section |4] deals with 
the analysis of the phase diagram of the binary collision model, by looking at the 
stationary homogeneous solutions and studying their linear stability. A transition 
toward collective motion is observed, but the spatially homogeneous motion turns 
out to be unstable in the validity domain of the hydrodynamic equations, namely 
close to the transition line. A comparison with the agent-based model is presented, 
showing that the transition lines of both models are qualitatively similar, and share 
some quantitative properties. Then, Section[5] investigates the behaviour of the binary 
collision model beyond the strict domain of validity of the hydrodynamic equations. 
A direct stability analysis shows that far from the transition line, the spatially 
homogeneous motion is stable. We then test whether the hydrodynamic equations 
could be used, in this domain, as a semi-quantitative description. We find that 
the restabilization phenomenon is indeed observed in the hydrodynamic equations, 
although the predicted location of the transition line between stable and unstable 
motion does not match quantitatively with a perturbative treatment of the kinetic 
theory. We also show that there exist solitary wave solutions of the hydrodynamic 
equations, that resemble the travelling stripes of higher density observed in the agent- 
based model. Finally, Section [5] discusses the main conclusions and perspectives of 
the present work. Some technical aspects related to the agent-based model and to 



Microscopic derivation of hydrodynamic equations for self-propelled particles 5 



the stability analysis of the homogeneous motion are reported in |Appendix A| and 
|Appondix B| respectively. 

3. Boltzmann approach and hydrodynamic equations 

3.1. Description in terms of Boltzmann equation 

One of the main goals of this work is to derive analytically from the microscopic 
dynamics, within an appropriate approximation scheme, the equations governing the 
evolution of the hydrodynamic fields, namely the density and velocity fields. A 
standard approach to obtain these hydrodynamic equations is to write, as a first 
step, the Boltzmann equation describing the evolution of the one-particle probability 
distribution in phase-space (i.e., the probability that a particle is at a given point, 
with a given velocity), and then to derive hydrodynamic equations by computing the 
first moments of the Boltzmann equation. Note however that such a procedure often 
yields a hierarchy of equations, so that a closure assumption has to be used. 

Let us start by deriving the Boltzmann equation for the above model. This 
equation relies on the standard assumption that the gas is diluted, meaning that the 
typical distance £ pp between particles is large compared to the interaction distance 
do, that is H 3> 1. In the present context, one also needs to assume that the ballistic 
distance ^bai is much larger than do, namely B ^> 1. It ensures that there is no 
memory effect from one collision to the other. The Boltzmann equation governs the 
evolution of the distribution f(r,9,t), that gives the probability that a particle is at 
point r with a velocity along the direction defined by the angle 9. On general grounds, 
this equation can be written as 
df 

^(r, 6, t) + v e(0) • V/(r, 6, t) = I dii [f] + I col [f, /]. (6) 

The different terms in the equation can be interpreted as follows. The second term in 
the l.h.s. corresponds to the ballistic motion of particles between two stochastic events 
(self-diffusion or collision). In the r.h.s., the term idif [/] accounts for the self-diffusion 
events, and it reads 

/7T />OC 
d9' / dnpoirj) (7) 
-7T J — OO 

OO 

x S(9' + n-6 + 2mTr)f{r,6',t). 

m— — oo 

The sum of ^-distributions accounts for the periodicity of angles. Finally, the term 
Icoi[f, f] describes the effects of collisions. It can be derived in the following way. A 
collision between two particles occurs if their relative distance becomes less than do- 
Although the two particles a priori play a symmetric role, it is convenient to choose 
one particle, and to observe the situation in the referential of this particle -say particle 
1. In this frame, the velocity of particle 2 is v 2 = v (e(9 2 ) — e(#i)). Hence, in order 
to collide with particle 1 between t and t + dt, particle 2 has to lie at time t (in the 
referential of particle 1) in a rectangle of length |v 2 |di and of width 2d . Coming 
back to the laboratory frame, this rectangle deforms into a parallelogram, but keeps 
the same surface, given by 2d o uo|e(0 2 ) — e(8i)\dt. The collision term I co i[f, f] is then 
obtained from the bilinear functional icoi[g, h]: 

Jcoi[<?,/i] = -2d v g(r,e,t) [ d6' |e(0') - e{6)\h(r, 6', t) (8) 
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+ 2d Q v d9 x \ d9 2 dr)p(v) |e(0 2 ) - e(0i)| 

J —TT J —TT J —OO 

OO 

X g(r,6i,t)h(r,02,t) ^ S(9 + ?] - 9 + 2rmr), 

Tfl — — OO 

with again the notation 9 = arg(e lSl + e 102 ) , and where g and h are arbitrary phase- 
space distributions. 

It is straightforward to check that the uniform onc-particlc distribution 
fo(r,9,t) — po/2n, associated to a uniform density of particles po, is a stationary 
solution of the Boltzmann equation, for any values of the noise parameters a 
and Co i since each term in Eq. ([6]) vanishes independently. If a transition to a 
state with collective motion occurs, another distribution should be a steady-state 
solution of the Boltzmann equation. Yet, finding this non-trivial distribution through 
non-perturbative analytical method is a hard task. One could turn to numerical 
approaches, but we would rather like to obtain analytical results, at least in some 
specific regime. We thus use in the following an alternative approach, which consists 
in deriving hydrodynamic equations for the density and velocity fields from the 
Boltzmann equation, in the limit of small hydrodynamic velocity. A stability analysis 
can then be performed on these hydrodynamic equations in order to check the onset 
of collective motion. 



3. 2. Derivation of the hydrodynamic equations 

3.2.1. Hydrodynamic fields and continuity equation. The hydrodynamic fields are on 
the one hand the density field: 

p{v,t) = ( d9f(v,9,t), (9) 



and on the other hand the velocity field: 

u(r, t) = -p-r- ^ d9 /(r, 9, t) e(9). (10) 

The equations governing the evolution of these hydrodynamic fields are derived by 
taking the successive moments of the Boltzmann equation. A simple integration of 
Eq. ([6]) over 9 directly leads the evolution equation for p{r,t): 

g+V.(pu)=0, (11) 

which is nothing but the usual continuity equation accounting for the conservation of 
the number of particles. 



3.2.2. Angular Fourier expansion of the phase-space distribution. The derivation of 
the evolution equation for the velocity field is actually much more complicated, and 
one has to resort to approximation schemes. As f(r,9,t) is a periodic function of 9, 
it is convenient to work with its Fourier series expansion, defined as: 

f k (r,t)= f d9f{v,9,t)e lke . (12) 
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Conversely, /(r, 0, t) can be expressed as a function of the Fourier coefficients through 
the relation: 

oo 

f(r,9,t) = - J2 Mr,t)e-* ke . (13) 

k— — oo 

In this framework, the uniform distribution /o(r, 0, t) = (2tt)~ 1 pq corresponds to 
/ fe (r,i) = (27r)-Vo4,o- 

Let us use as a basis of the plane the two orthogonal vectors ei = e(0) and 
e2 = e(n/2). The components of e(0) in this basis are obviously e\(9) — cos0 and 
62(0) — sin0. In order to obtain an evolution equation for the velocity field, we 
multiply Eq. © by e(0) and integrate over 9; one gets in tensorial notations (j = 1, 
2): 

^ J" ie^Oj/fr.Ml+^^j 1 d9e ] (9)e l (9)f(r,9,t) = (14) 

^(0) (7 di f[/]+ /col [/,/]). 

To proceed further, it is convenient to identify complex numbers with two-dimensional 
vectors, in such a way that e(9) is mapped onto e 10 . Then, in the same way, i>o/i(r, t) 
is associated to the momentum field w(r, t) — p(r, t) u(r, t). Hence, we wish to rewrite 
Eq. (fT4)l in such complex notations. For later use, we shall write it in a slightly more 
general form, replacing e lB with e tke (k being an integer): 

2 

d9e lke e l (9)f{v,9,t) = (15) 

1/. 1 1 i _ 



f ^e lfee /(r,0^)+ V oE7T- 1 

./— 7T ./— 7T 



d0e ifc " (Idif[/] + /coi[/, /])• 



Eq. I|14p is recovered for fc = 1, up to the mapping between complex number and 
two-dimensional vectors. The first term in the l.h.s. is simply dfk/dt. The r.h.s. of 
Eq. (fT5|) is computed by inserting the Fourier series expansion Q13[) into Eqs. ([7]) and 
(|5J). After a rather straightforward calculation, one finds: 



dO 



ikO 



(JdifLf] +Icol[/]) = -A (l - e^' 2 ) / fc ( r>t ) (16) 



2 



-d Q v Y, (l q -er k ° l 2 I q _ k/2 ) f q {v,t)f k ^ q {v,t) 



q— — oo 



where the coefficients I q are defined as: 



h = 





. 9 


f d9 


sm — 


I — 7T 


2 



cosg0. (17) 



From this definition, it is obvious that I- q = I q . For integer q, I q is given by: 

= I^V ' (18) 

while for half-integer g = m + 5 (m integer) one has: 

/i=I_ | =2, (19) 

Z m+ x =_L_[(-l)'»(2m+l)-l] 1 m^-1,0. (20) 
2 m(m + 1) 
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The second term in the l.h.s. of Eq. (To)) can be evaluated as follows. For 1 — 1,2 and 
k integer, let us define the complex quantity Q; (r, t) as: 

Q { ?\v,t)= r dee ike ei (e)f(r,9,t). (21) 

J — TT 

The following relations are then easily obtained: 

Q[ k) (r,t) = ~[/ fe+ i(r,t) + /*_!(!■,*)], (22) 

Q?\r,t) = l[f k+1 (r,t)- f k -x(r,t)}. (23) 

3.2.3. Velocity field equation in the small velocity regime. Up to now, the calculations 
made are exact, apart from the approximations underlying the Boltzmann equation. 
As already mentioned, the Fourier coefficient /o(r, t) is nothing but the density field 
p(r,t), and /i(r, t) can be mapped onto the momentum field w(r,t) through the 
identification of complex numbers with two-dimensional vectors. A similar mapping 
also holds for /_i(r, t), which is the complex conjugate of /i(r, t). In contrast, Fourier 
coefficient /^(r, t) with |fc| > 1 cannot be mapped onto the hydrodynamic fields. As it 
turns out that such coefficients appear both in the expression of <3, (r, t) and in the 
r.h.s. of Eq. (|15p . an approximation scheme has to be found in order to obtain from 
Eq. (TK)|) a closed hydrodynamic equation, involving only the fields p(r,t) and u(r,i). 

In the following, we assume that the distribution /(r, 9, t) is close to an isotropic 
distribution, namely, it depends only slightly on 9. This amounts to assuming that 
the hydrodynamic velocity is much smaller than the velocity of individual particles. 
In terms of Fourier coefficients, the hydrodynamic velocity is given by ||u(r, t)\\ — 
i>o|/i(r, i)|/p(r, t). We introduce a small parameter e such that ||u(r,t)|| = 0(e). For 
instance, e can be chosen as u/vq, where u is the spatial average of ||u(r, i)|| at some 
initial time t = tg. Then the key assumption we use to build an approximation scheme 
is 

A(r,i) -C(el fc l). (24) 

Such a scaling ansatz is consistent with the property /_fc(r,i) — fk(r,t)*, with the 
scaling properties of /o(r,i) and /i(r,i), and with Eq. (fTBjl . We shall identify more 
precisely in Section 14.1.21 the validity domain of this scaling ansatz, and thus of the 
hydrodynamic equations we will derive from it. 

Using the above scaling ansatz, the sum in the r.h.s. of Eq. (jT6j) , for k = 1, can 
be truncated, only keeping terms with q = 0, 1 or 2, that are at most of order e 3 , 
while discarding the other terms, being of higher order in e. Gathering all terms, one 
obtains the following equation for the evolution of f\ (we drop the explicit dependence 
upon r and t to simplify the notations): 

dfi v d * v d f 



A (l - er^) + V g - e- 2 / 2 ) 



-dovo ( e-° 2 / 2 - \ ) flh- (25) 



7T V 5 
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Hence, the resulting equation involves /o = p, f\ and / 2 . Accordingly, it turns out 
that one needs to find a closure relation to express fi as a function of /o and fx (or, 
equivalently, in terms of p and u) . Such a relation is given by the evolution equation 
for /2, that is, Eq. (|T5)) with k = 2. From Eq. dTHJ) , one sees that Fourier coefficients f q 
with |g| > 2 are a priori involved, but they can actually be discarded as being of order 
higher than e 2 , whereas f 2 = 0{e 2 ). Similarly, the quantity Q^p can be expressed as 
a function of f\ and /3, and here again, f-$ can be neglected. One thus ends up with 
the following equation for fi'- 



dt 



vodfi 
2 dx! 



vq dfi 
2i dx2 

A (l - e 



-dov - 



16 A 
6TT 



ft- 



(26) 



Within our hydrodynamic description, it is also natural to assume that the phase- 
space probability density /(r, 9, t), or equivalently, its Fourier coefficients /fc(r, t), vary 
significantly only over time and length scales that are much larger than the microscopic 
ones. Relevant microscopic time scales are the typical collision time r co i = l/(pdova), 
and the typical ballistic time Tbai = 1/A between self-diffusion events. It is thus 
legitimate to neglect the term dfijdt in Eq. (|26|) . as it is much smaller than fijT co \ 
and /2/Tbai- In contrast, the terms containing the spatial derivatives have to be 
retained, as they involve f\ which is much larger than fi. 

From Eq. (|26|) -without the time-derivative term- one can express fi as a function 
of p and fx. Then plugging this expression for fi into Eq. (|15|) . with k = 1, leads 
to a closed hydrodynamic equation governing the evolution of fx, and involving only 
fx and p. Mapping back complex numbers onto two-dimensional vectors, vofx can 
be identified with the "momentum" field w = pu, and one obtains the following 
hydrodynamic equation: 



— + 7(w • V)w = - yVp + - Vw 2 + (p - £w 2 )w - 



iA7 2 w 



k(V-w)w + 2zA7 / o-M-i/(V-w)V / 9, (27) 
Vw T ) is the symmetric part of the 



with v' = dv/dp, and where M = |(Vw 
momentum gradient tensor. The different coefficients appearing in this equation are 
given by: 



4 

16^do 

TTVq 

l6isdo 

TTVQ 

8 . 

-d v p 
ir 



A 1 1 — e" 



— + 2e 



2e" 



-a 2 /2 



-a 2 /2 



16 A 

i^-d Q v p 
Sir 



7 



,-o>/2 

■X(l 

1 

3 



-ag/2 



-2(7' 



(28) 
(29) 
(30) 

(31) 
(32) 
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Eq. I|27p may be considered as a generalization of the Navier-Stokes equation to a 
case where on the one hand, the global momentum of the assembly of particles is 
not conserved by the microscopic dynamics, and on the other hand, the dynamics 
breaks the Galilean invariance. This shows up in the appearance of new terms in the 
equation, as well as in the presence of the coefficient 7, generically different from the 
Navier-Stokes value 1/p, in front of the (w • V)w term. For instance, if A <C pdovo, 
7P remains close to 0.6 for any value of a. 

The different terms in the r.h.s. of Eq. (|27|) may be interpreted as follows. 
Neglecting the density dependence of k, the first two terms can be considered as 
a pressure gradient, where the effective pressure P e e obeys the equation of state 
Pes = \{vqP~ kw 2 ). The third term accounts for the local relaxation of the momentum 
field w, and this term plays an important role in the onset of a collective behaviour, 
as we shall see in the following section (note that £ > when p > 0). The fourth 
term describes the viscous damping, like in the usual Navier-Stokes equation. The 
parameter v can thus be interpreted as a kinematic viscosity. It decreases when p 
increases, but the 'dynamic' viscosity pv increases with p. The fifth term may be 
thought of as a nonlinear feedback on the momentum field of the compressibility of 
the flow. Finally, the two last terms correspond to a coupling between the density and 
momentum gradients. 

It is also important to note that the above hydrodynamic equation (|27[) is 
consistent with the phenomenological equation postulated by Toner and Tu on the 
basis of symmetry considerations 17:. Specifically, expanding the expression of 
w = pu in that equation, we find the same terms involving the velocity gradients 
as in Ref. [17]. But it turns out that the term V(V • w), that would be allowed 
from symmetry considerations, does not appear in the present approach, that is, the 
coefficient in front of it vanishes. Note also that the term (u • V) 2 u considered by 
Toner and Tu [17j . does not appear here for being of higher order than the terms 
retained in the expansion. Last, hydrodynamic equations which have been derived 
through the kinetic approach are entirely deterministic, while Toner and Tu studied 
stochastic equations. However some additional terms also appear, like the coupling 
terms between density and velocity gradients. Most importantly, the present approach 
provides a microscopic justification to the hydrodynamic equation of motion, and 
yields explicit expressions, as a function of the microscopic parameters, for the different 
coefficients appearing in the equations (transport coefficients). 

4. Noise-density phase diagram from the hydrodynamic equations 

Spatially homogeneous stationary solutions 

4-1.1. Transition toward collective motion. Now that the hydrodynamic equations 
of motion have been derived, it is natural to look for the different possible stationary 
solutions and to test their stability. Let us first look for the spatially homogeneous 
stationary solutions. Dropping all space and time derivatives, one ends up with the 
simple equation: 

(/j - £w 2 )w = 0. (33) 

Hence a trivial homogeneous stationary solution is w = for all values of the 
parameters. When p > 0, a second solution appears, namely w = Wi = \fp~R, e ; 
where e is a unit vector pointing in an arbitrary direction. The stability against 
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Figure 1. (a) Phase diagram of the model in the plane (p,cr), with p = pvodo/X. 
A transition line (full line: <tq = <r; dashed line: erg = 1) indicates the linear 
instability threshold of the state u = |u| = 0. (b) Hydrodynamic velocity u in the 
homogeneous state for a = <ro = 0.6, computed numerically from the Boltzmann 
equation (full line) and analytically from the hydrodynamic equations (dashed 
line). Inset: same data on logarithmic scales (dots: slope 1/2). 




spatially homogeneous perturbations is easily tested by assuming that the flow is 
homogeneous, but time-dependent in Eq. (|7f|) , yielding: 

— = ( M - £w 2 )w. (34) 

It follows that w = is a stable solution when p, < 0, while it becomes unstable for 
p > 0. In the latter case, the emerging solution w = Wj is stable against homogeneous 
perturbations. From the expression ([3"T]) of /i, we see that the sign of p is related to 
a competition between density and self-diffusion. When the self-diffusion probability 
A is high, p < and there is no flow. In constrast, when the density is high, p > 
and a spontaneous flow appears, due to the numerous interactions between particles. 
The value = defines a transition line in the phase diagram noise versus density: 
for given values a and do of the noises, the nonzero solution w = wi appears for a 
density p > p t , where the threshold density p t is given by: 

ttA(1 - e- ff o/2) 

8d wo(e ' - f) 

In terms of the dimensionless parameter (or reduced density) 

B d 4ai pdovo 
P= JP=^ = — (36) 

the threshold is expressed as 

Pt = 8(e— (37) 

This last result is interesting, as it shows that the threshold pt, which could a priori 
depend on the three dimensionless numbers er, do and B, actually does not depend 
on B. The transition line is plotted in Fig. []Ja) for the two cases o~o = a and <To = 1. 
Instead of considering the transition as a function of the density, one may also look for 
the transition by varying the noises at a given fixed density. If the two noise intensities 
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do and a are equal, the instability of w = occurs for any (non-zero) density, and the 
threshold noise cr t behaves in the low density limit p — > as <r t ~ p 1 ! 2 . This nontrivial 
prediction can be verified in direct numerical simulations (see below). In contrast, 
when (To is kept fixed while varying <r, no transition occurs as a function of a if the 
reduced density is lower than a limit p® given by: 

P° t =^(l-e-^). (38) 

Finally, in the opposite limit of high density, the threshold noise er t saturates to a 
value of = (2 In f) 1 / 2 « 0.90. 

4-1-2. Validity domain of the hydrodynamic equations. The hydrodynamic equations 
rely on the scaling ansatz (|24|) . In order to verify a posteriori the validity of the 
hydrodynamic equations, we compare the stationary homogeneous solutions with 
non-zero velocity obtained from the hydrodynamic equations to that numerically 
computed from the Boltzmann equation. The hydrodynamic velocity u, computed 
as ti = u\ = p^ 1 \J /V£> i s plotted on Fig. [IJb) as a function of the reduced density 
p. Note that u/vq is a function of the dimensionless numbers p, a and <7o only. As 
expected, the velocity u computed from the hydrodynamic equation matches perfectly, 
in the small velocity regime (i.e., close to the transition line) the numerical data from 
the Boltzmann equation. However, it turns out that even quite far from the transition, 
when u becomes of the order of vo, the value u\ computed from the hydrodynamic 
equation remains a good estimate of the value obtained from the Boltzmann equation. 
In particular, it is interesting to note that u\ also saturates at large p to a finite value 
uf(a) < 1, given by 



u-{a)=v 



' 1 -\ (ti + - 



3 M 15 1 3 



(39) 



(^ 2/2 -I)(l + ^ 2 ) 

(see Fig. Gib)). Hence, even beyond their domain of validity, which is restricted 
to small values of the hydrodynamic velocity, the hydrodynamic equations we have 
derived yield a rather good approximation of the exact dynamics. Specifically, they 
fulfill the condition that the hydrodynamic velocity should remain smaller than the 
individual velocity vq of the particles, although this result was not a priori obvious 
given the approximations made. 

To further test the validity of the hydrodynamic equations, we have also checked 
explicitely, from a numerical calculation, that the scaling ansatz (|24l) is correct. 
Specifically, we computed from a numerical integration the stationary and spatially 
homogeneous solution of the Boltzmann equation. In order to work with 
dimensionless quantities, we plot on Fig. [U(a) the quantities gu = /jf// 3 (instead of 
/I*) as a function of k. We observe that gu decays almost exponentially with k, as 
soon as k > 4. To test the scaling ansatz, we first reformulate it in a more specific 
way. The ansatz is obeyed if there exists for all k a constant c& such that gu ~ Ckg\ 
in a parameter regime where g\ -C 1. We thus plot on Fig. [U(b) the ratio gkj g\ for 
different values of the density, close to the transition, and we observe a reasonable 
collapse of the data. Let us however emphasize that a strict collapse of the data is 
not necessary in order to apply the approximation scheme used in the derivation of 
the hydrodynamic equations. The essential requirement is that the quantities gk with 
k > 2 could be neglected. As the ratio gkjg\ decays rapidly with k, neglecting terms 
with k > 2 is a safe approximation. 
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Figure 2. Test of the scaling ansatz /j. = ). (a) = ftf/p versus k, for 

cr = cro = 0.6 and different values of the reduced density p, close to the transition 
density (pt = 0.3837); an exponential decay is observed, (b) g^ldX as a function 
of k, showing that for a given k, is essentially proportional to g£ when the 
density is varied (A = 0.5, dg = 0.5, vq = 1). Inset: zoom on the small k region. 



4-2. Stability against inhomogeneous perturbations of the homogeneous stationary 
solutions 

4-2.1. Evolution equation for the perturbations. We have shown that above a 
threshold density p t , or equivalcntly, below a threshold noise at, the solution with zero 
velocity becomes unstable, and a stable solution with finite velocity emerges. Yet, only 
the stability with respect to homogeneous perturbations (i.e., with infinite wavelength) 
has been tested up to now. Hence this does not ensure that the finite velocity solution 
is really stable, as it may be destabilized by finite wavelength perturbations. We 
now check this issue, by introducing small perturbations around the homogeneous 
stationary solutions po and Wo, namely 

P (r,t) = Po + 8p(r,t), w(r,t) = w + <5w(r, i). (40) 

Note that Wn may either be equal to zero or to the nonzero solution wi. Plugging 
these expressions into the hydrodynamic equations (jlip and (|27|) , we can expand the 
resulting equations to first order in the perturbation fields 6p(r,t) and <5w(r, t), also 
taking into account the density dependence of the different coefficients. This yields 
the following linearized equations: 

d 

— 5p + V-6w = Q. (41) 
at 

d v 2 

— Sw + 7(w • V)Sw = — + «V(w • Sw) (42) 

+ [(// — £'wg)(5/9 — 2£wo • Sw — kV • 5w]wo 
+ (ji — £wg)<5w + ^\7 2 <5w, 

where p! and £' are shorthand notations for dp/ dp and d^/dp. Note that dp/ dp is 
actually a constant, i.e., it is independent of p. Then we make the following ansatz 

5p(r , t) = Sp e st+lcl - r , Sw(r , t) = <5w e st+l * r , (43) 

where Smvq is a vector (with real components), and Spo is a complex amplitude that 
takes into account a possible phase shift between density and momentum perturbation 
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fields. Both ||<5w || and \5po\ are assumed to be small. The wavenumber q is assumed 
to have real components, whereas the growth rate s is a priori complex. In addition, 
q is considered to be given, and one looks for the dispersion relation s(q). If the real 
part Sft[s(q)] > 0, the mode with wavenumber q is unstable. Then Eqs. (|TTj) and ([4*2")) 
become: 

s Sp + iq ■ <5w = 0, (44) 
[s + 7(w • iq) - (p - £w„) + ^q 2 ]c$w = 

- 2^0 S Po ~ 2kw o • <5w )iq 

+ [{p! - C'wg)(5p - (2£w + niq) ■ 5w ]w . (45) 

Note that, due to linearity, the above equations can be re-expressed as a function of 
the ratio of amplitudes 8w Q /SpQ. 



4-2.2. Stability of the zero-velocity solution. Let us first check the stability against 
inhomogeneous perturbations of the solution Wo = 0, which is known to be stable 
against homogeneous perturbations in the low density phase p < pt (corresponding to 
p < 0). In this case, Eq. (|45jl simplifies to: 

(s + vq 2 - p)Sw = --v$5p q, (46) 

Thus <5wo is along the same direction as q. Writing q — qe and <5wo = Swoe, where e 
is an arbitrary unit vector, one can eliminate the ratio Swq /Spo from Eq. (|44p , yielding: 

.sS + OV-^ + lY = 0. (47) 

The discriminant of this second order polynomial equation reads (note that p < 0): 

A = (\p\+isq 2 ) 2 -2v 2 q 2 . (48) 
If A > 0, the roots are real, and one finds for the largest one s + : 
1 



-(\p\ + vq 2 ) + ^/(\p\ + vq 2 ) 2 - 2v 2 q 2 



< 0. (49) 



In the opposite case A < 0, the roots s± are complex conjugates, and their real part 
is given by: 

K[s±] = -i(|Ml+^ 2 )<0. (50) 

As a consequence, the homogeneous fields wo = is stable with respect to finite 
wavelength perturbations in the region p < p t . 



4-2.3. Stability of homogeneous collective motion. We now turn to the stability 
analysis of the stationary homogeneous flow Wo = wi, obtained for p > p t . For 
the hydrodynamic equations to be valid, we restrict our study to values of p very close 
to pt, with p > p t . One could a priori consider vectors q and 5\vq that make arbitrary 
angles with respect to Wi. However, it can be shown (see |Appcndix B| that only 
some specific angles are allowed. Further, for all allowed perturbation modes such 
that q and <5wo are not along the direction of wi, the real part of the growth rate s 
is negative, so that these modes are stable ( |Appcndix B[ ). The only instability that 
appears is for longitudinal perturbations, such that q, Swq and Wi all have the same 
direction. We thus focus on this specific case in the following. 
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Figure 3. Longitudinal instability, (a) 5R[s_|_]/A versus q for p = 0.22 (full line), 
0.30 (dashed) and p = 0.4 (dot-dashed), and a = a = 0.5 (p t = 0.2138). The 
maximum growth rate decreases when p is increased, (b) (qjBdo) vs (p — pt) in 
logarithmic scales. The dashed line indicates the scaling g,; <x (p — ■ pt) 1//4 . 



Considering a longitudinal perturbation, we write Wi = w±e, q = qe and 
Swo = Swoe, where e is a unit vector. Under these assumptions, Eqs. (j44| and 
Toll become: 



s Spo + iq Swq = 0, 
(s + ■yiqwx + vq 2 )8wo 

+ Wi[(fl'- 



"1 



(51) 
(52) 



£'wl)6p - (2£wi + iqK)6w a ], 



where we also take into account that p — = 0. From Eq. ([51]) , one gets 



Swo/dpo = —s/iq, which we report in Eq. ((5 
degree in s 

„2 , . r/.._2 



[(vq 2 + 2p) + iq^ywi] 



r 2 2 

r ?v 



+ iqwi (p - gwl) 



This yields a polynomial of second 

(53) 

= 0, 



from which two solutions s± can be obtained. Denoting as s + the solution with the 
largest real part, we find 



K[ S+ ] = -~K 




(54) 



with 



Ji 



{vq 2 + 2p) - q 2 (^w 2 + 2v 2 ) 
2 Wl q[-f(vq 2 + 2p) - 2p' + 



To deal with this complicated expression, we first plot 3?[s+] as a function of q for some 
specific values of the parameters (see Fig. G^a)). Near the threshold p t of collective 
motion, there exists a threshold value qi such that 5ft [s+] is positive for q < qi and 
negative for q > qi. Hence the homegeneous flow turns out to be unstable with respect 
to long wavelength perturbations. 
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This result is confirmed by a small q expansion of Eq. (|54p . Expanding 9?[s+] up 
to second non-trivial order in q, that is to order q since only even powers of iq appear 
in the expansion of the real part of s + , we gefl 

The positivity of the coefficient of the q 2 term confirms that, close to the transition 
line, long wavelength modes are unstable. Note that the expansion (|55|) is consistent 
as long as the fourth order term remains small with respect to the second order one, 
yielding the condition q <C q* which defines the wavenumber q*: 

q* =d^ 1 B- 1 (p- Pt )H(p). (56) 

The function ^(p) goes to a constant value for p — > 0, and ~ p^ 1 / 2 for p — > oo. 
The wave vector q* defines the region where the first term of the expansion of 5i[s+] 
is dominant. 

Interestingly, we observe that other wavenumbers characterizing 9?[s+] have a 
different scaling with p — pt, the deviation from the threshold. For instance, it 
can be shown analytically that the wavenumber qi, defined by 3?[s+] = 0, scales as 
qi ~ w\^ 2 ~ (p — Pt) 1 ^, as illustrated on Fig. [2Kb). The wavenumber delimitates 
the domain of unstable modes. Another example is given by the wavenumber q m that 
maximizes 3?[s+], and thus corresponds to the most unstable modes: q m is found to 
scale as q m ~ W\ ~ (p — pt) 1 ^ 2 - The existence of these different scaling regimes is an 
illustration of the complexity of the dynamics close to the transition line. 

Finally, we emphasize that the perturbations that destabilize the homogeneous 
collective motion (that is, the long-range order) are different from the ones that 
destabilize long-range order in the XY-model, an equilibrium model with essentially 
the same symmetries as in the present model. In our model, motion is destabilized 
by longitudinal waves, while in the XY-model, long-range order is destabilized by 
spin-waves, that is, by a small change in the spin direction from one spin to the 
neighbouring ones. 



4-3. Comparison with the phase diagram of the agent-based model 



4-3.1. Numerics and parameters. All simulations are performed using models defined 
on a square domain, with periodic conditions on both boundaries. The initial 
conditions always consist in randomly dispersed particles, with a uniformly chosen 
random speed direction. Then all measurements are performed after a sufficiently 
long time so that a stationary state is reached. 

In the above framework of the Boltzmann equation, we considered diluted 
systems with small correlations between particles, which is expressed in terms of the 
dimensionless number H and B as: 



H > 1, 



B> 1. 



(57) 



X To simplify the resulting expressions, we approximate the coefficients of the expansion in q by their 
leading order in l//t, as fi is small close to the transition line. The full expression of the coefficient 
S2 of the q 2 term reads: 



•7 



2_ 



2vl 



This expression will be used in Fig.[5]to compare with numerical results. 
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In the numerical agent-based model, we do not have access to very large values of H 
and B, due to simulation constraints. However, to be as consistent as possible with the 
kinetic theory approach, we mainly explored a parameter range such that H > 4 and 
B > 4. Among the three dependent dimensionless numbers B, H and p, we decided 
to keep B to characterize the set of parameters, and p as the control parameter. 

Throughout the study, we fix At — 1. We defined some sets of parameters 
(do,vo,X) and, for each of them, we studied the behaviour of the system in the 
parameters space (p,r]). To make the comparison between analytical and numerical 
results easier, we characterize the noise amplitude by its rms-value a (or equivalently 
its variance a 2 ). For a uniform noise on the interval [— 7771", 7771"], we have a = r\.^j\fi. 
The self-diffusion noise is kept equal to the collision noise (cto = a), except for one set 
of parameters in which the angle of diffusion 770 £ \ is chosen over the whole circle. We 
call its rms-value cr™ ax . All the parameter values are summarised in Table [T] 

4-3.2. Transition line. When the noise amplitudes for collision and self-diffusion 
are equal, the general aspect of the phase diagram is the same both for the kinetic 
theory, and the agent-based model (Fig. Ufa)). We have drawn the transition line for 
different sets of parameters on Figure 0Jb) . All curves seem to be bounded between 
configurations I (B = 0.5) and II (B = 4). 

Looking at the influence of the different parameters, we can make the following 
observations. First, there are small variations as the self-diffusion probability A 
changes with a factor of eight from configuration II to IV (Fig. BJb)). When the 
dimensionless parameter B is kept constant (Sets V and VI), the transition points 
corresponding to the same value of p are equal within the error bars. Apart from Set 
I for which B < 1, it turns out that the measured values of at differ by less than 15% 
for any given p, while B is varied by a factor of 16 between Set II (B = 4) and Set VII 
(B = 64). However, we are not able to conclude that the curves merge into a single 
master curve. In particular, the observed evolution of at when increasing B at fixed 
p is not monotonous (Fig. BJb)). 

In the low p region, the transition noise varies as a power law with p, a oc p@ when 
p — > 0. We have measured the exponent (3 for the largest dimensionless number B 
(B = 64, set VII), yielding (3 = 0.46 ± 0.04 (Fig. @fc)). This value is compatible with 
a square-root behaviour as found analytically in the binary collision model (Fig. [T]). 
Quantitatively, the transition line computed from the kinetic approach and the one 
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Table 1. Physical parameters of agent-based simulations; set I corresponds to 
the values used in [15] . Parameters are chosen as multiple or sub-multiple of 2; 
g-max _ Tr/y^ corresponds to the same variance as a uniform noise on [— it; it]. 
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Figure 4. Phase diagram of agent-based models, (a) Overview of the phase 
diagram. The continuous line is the transition line of the continuous model given 
in Eq. H37I I; data for symbols o are obtained with the parameter set I. (b) Diagrams 
for all configurations with the self-diffusion noise co = a (plot in log-log scales), 
(c) Scaling of the transition line at small p with the set of parameters VII (B = 64) 
in log-log scales. The continuous line is the transition line H37H obtained from the 
kinetic theory. The dashed line is a square-root fit of the numerical results, (d) 
Model with a constant and maximum noise amplitude for self-diffusion in log-log 
scales (set VIII, B = 64). The continuous line corresponds to a fit of the numerical 
points with the law at = a(p — pj 1 ) 1 / 2 . The dashed line is a fit with a power law 
/3 = 0.57 (see Table[T]for values of the other parameters). 

which we measure in the agent-based model are relatively close one to the other; their 
largest relative difference is about 30%. 

4-3.3. Maximal self- diffusion. When we set the amplitude of self-diffusion noise to 
its maximum (770 = 1 or ctq = w/yS), the behaviour of the model remains qualitatively 
similar to the case a = o~q that we studied above, with only a few quantitative 
differences. The transition line is shifted to a lower noise amplitude: ct differs by two 
orders of magnitude between the two comparable parameter sets VII and VIII. Fitting 
the two curves by a power law, the exponents are significantly different: [3 ~ 0.46 for 
set VII (cr = er ), while j3 « 0.57 for set VIII (cto = crJJ lax ). One possible explanation 
for such a difference would be that, as in the hydrodynamic equations, there exists a 
threshold pj? below which no collective motion occurs, whatever the noise amplitude 
a. A fit with the function p = a\Jp — p° gives a value p° = 0.00133 (Fig. HJd)), much 
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smaller than the theoretical value — 3n/8 w 1.18. Given the presently available 
data, we are not able to discriminate between the two fits, and to conclude on the 
existence of a non-zero threshold value pj?. Trying to find a phase transition for a 
very low value of p (p = 2~ 10 ), below the fitted value p®, we could hardly define a 
threshold. However, it might be necessary to reach larger system sizes to detect a 
phase transition in this regime. 



5. Beyond the strict validity domain of the hydrodynamic equations 

In Section [4l we concluded from a linear stability analysis that the homogeneous flow 
is unstable with respect to long wavelength perturbations, in the validity domain 
of the hydrodynamic equations, namely close to the transition line. When getting 
farther from the transition line, previous theoretical approaches |18] suggest that 
the homogeneous motion should be stable. To come to a conclusion in the present 
framework, it is thus necessary to come back to an analysis of the Boltzmann equation. 
It is also natural to wonder whether the hydrodynamic equations could yield, out of 
their strict validity domain, a qualitative description of the phenomenology of the 
moving phase. We address these issues in the present section. We find in particular a 
restabilization of the homogeneous flow far from the transition line, as well as solitary 
waves that we compare with the travelling stripes already reported in numerical 
simulations of the agent-based model [14] . 



5.1. Stability analysis from the Boltzmann equation 

In order to analyse the stability of the finite velocity solution beyond the validity 
domain of the hydrodynamic equations, we come back to the Boltzmann equation, 
and we resort to a semi-analytical treatment. 

We start with a formal expansion of the phase-space distribution /(r, 9, t) around 
the homogeneous stationary solution fo(9): 

f(v,e,t) = fo(d) + 6f(r,e,t). (58) 

Considering a perturbation of wavevector q of the form: 

6f(r,e,t) =<5poG(0,q)e s * +iq - r , (59) 

with d0G(6,q) = 1. Assuming, as in Sect. 14.2.31 that both q and the velocity 
perturbation are along the same direction e as the collective velocity, the function 
G(9, q) satisfies the following linearized Boltzmann equation: 

sG(0,q) + iqv cos 9 G(9,q) = I dii [G] + I col [G, f ] + / col [/ , G}. (60) 

Setting q = ge, we are interested in a small q expansion of Eq. (|60p . in order to 
compare with the results of Eq. (|55|) . We then expand s and G(8,q) in the following 
way: 

s = i Sl q + s 2 q 2 +<D(q 3 ), (61) 
G(0,q) ^G (9)+tqG 1 (9)+q 2 G 2 (9)+O(q 3 ), (62) 
with the normalization conditions: 

/TT pTT pTT 

d9G (9) = l, / d9G 1 (9) = 0, / d9G 2 (9) = 0. (63) 
-7T J —TT J —TT 
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Then Go, G\ and Gi are solutions of the hierarchy of equations: 

/dif [Go] + led [Go, f ] + /coi[/o, Go] = 0, (64) 

Sl G {9) + v cos6G (6) = 7 dif [Gi] + / co i[Gi, / ] + I C oi[/o, Gi], (65) 

-«iGi(0) + a a Go(0)-«OCOsflGi(0) = (66) 

^dif [Ga] + / C ol[G2, /o] + X=ol[/o> G2]. 

Using the properties of /dif and i co i, namely 



dBIm\g]=0, / dei col [gj ]= dei col [f ,g] = (67) 
for any function <?, we obtain 

Si = -u / d(9 cos6»G (6i), (68) 



s 2 = uo / d(9 cos0Gi(0). (69) 

Hence the determination of G2 is not necessary to compute S2- We only need to 
compute the hierarchy of functions up to G\. It is actually convenient to work in 
Fourier space, introducing the Fourier series expansion Gg,k and G\^ of Gq(9) and 
G\{0) respectively. In terms of this Fourier expansion, one finds s% = — Go,fc=i and 
s 2 = G\ t k=\i assuming that Gq{6) and Gi(6) are even functions. 

The integral equations (|64|) and (|65|) can be solved numerically, once expressed in 
terms of Fourier coefficients. To this purpose, we use the following Fourier expansion 
of the integral operators /dif and J co i: 

d6e m I An \g] = - A (l- e^**' 2 ) g k , (70) 
r d ee^I col [g,h} = 2 -^ £ - I q ) g k -X- (71) 

ff q— — 00 

Numerical results are reported in Fig. [5ja), where S2 is shown as a function of a for 
a = (To j all other parameters being kept fixed. Consistently with the results obtained 
from the hydrodynamic equations, we observe that close to the transition line, S2 is 
positive and diverging. But for smaller values of the noise amplitude <r, S2 becomes 
negative. Hence the homogeneous state of motion becomes stable in this range with 
respect to long wavelength perturbations. 

As summarized on Fig. EJb), there are from the point of view of stability three 
regions in the phase diagram (we focus here on the case a = <7o)- These three regions 
can be described as follows: 

A At low p or high a, no collective motion occurs. 

B For pt < p < p r , a homogeneous stationary solution with nonzero velocity exits, 
but it is unstable under longitudinal compression modes. 

C For p > p T , the homogeneous and stationary moving phase is linearly stable under 
any small perturbation. 

p T is defined as the value of the reduced density such that S2 — 0. Note that p T 
is not a monotonous function of a. In the B region, the system cannot converge 
to a homogeneous stationary solution, and one thus expects the system to organize 
into more complicated spatio-temporal structures, that we shall try to describe in 
Section 15.31 
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Figure 5. (a) Dependence of S2 on a for p = 2.5 and cto = u. For a below a 
given threshold cr r , S2 becomes negative, indicating that the homogeneous state of 
motion is stable with respect to long wavelength perturbations. In contrast, close 
to the transition line, this state is unstable since S2 > 0. The vertical dashed 
line corresponds to the transition value at- Inset: comparison, close to at, of 
S2 obtained numerically from the Boltzmann equation (full line) and analytically 
from the hydrodynamic equations (dashed line), showing a good agreement, (b) 
Phase diagram indicating, for a = erg, the three different regions: no motion 
(A), unstable homogeneous motion (B), stable homogeneous motion (C). The 
full line has been obtained numerically from a stability analysis of the Boltzmann 
equation. The dashed one is the transition line shown in Fig. \V[ a). 



5.2. Restabilization of the homogeneous flow in the hydrodynamic equations 

The above stability analysis from the Boltzmann equation shows that the homogeneous 
flow becomes linearly stable when getting farther from the transition line p t . Although 
this region of restabilization is, strictly speaking, out of the validity domain of the 
hydrodynamic equations, it would be interesting to know whether these equations 
already contain, at a qualitative level of description, the restabilization phenomenon. 

One possible way to investigate this stability issue is study the sign of the 
coefficient S2 of the q 2 term in the small q expansion of 5R[s + ] . An equivalent procedure, 
that we follow here, is to look for the domain of existence of the wavenumber qi (defined 
as 5R[s+] = for q. L ^ 0), when the control parameter p is increased at a given noise 
amplitude a. In order to achieve this task, we solve the equation 5R[s+] = 0, using 
expression ()54|) . The solutions are naturally expressed in terms of the variable qf. 
After some algebra, we find for the largest solution: 



2 A 4 
Qi = —2 



where the term y^- — ^-J is positive. The expression in the right hand side of Eq. (|72|) 

is positive for p close enough to p t , and becomes negative for larger p (see Fig. [(^a)), 
in which case a real solution qi does not exists. As a result, there exists a value p r of 
the control parameter p such that qi is no longer defined. For p > p Xl 9?[s+] remains 
negative for all values of q (Fig. [3ja)), so that all perturbations are linearly stable. 
Using equations (I28p -(|32 |) . we can compute the restabilization line p r (c) and show 
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Figure 6. Rcstabilization in the hydrodynamic framework, (a) (qiBdo) 2 such 
that 5R[s] = versus p, same parameters as Fig. \3\ (b) Phase diagram. The full 
line corresponds to the onset of motion, at- The dashed line is the transition 
between stable and unstable homogeneous flows, cr r . Regions A: w = 0, B: w ^ 
and Sft[s] > when q < qj, C: w ^ and R[a] < for all q and for all direction of 
perturbation. Inset: same as (b) in log-log scale. 



that p r depends only on a and oo; but not on B. We also find that p T (cr) behaves for 
small noise amplitude as p T oc c 1 / 2 (see Fig. Ufb) and its inset). 

Altogether, the hydrodynamic equations seem to lead to the correct 
phenomenology even when used beyond their strict validity domain. Yet, the locations 
of the transition line p r (o~) predicted from the hydrodynamic equations on one side, and 
the one predicted from a long wavelength perturbative treatment of the Boltzmann 
equation on the other side are quantitatively different, as illustrated on Figs. Efb) 
and^b). 

5.3. Inhomogeneous travelling solutions and solitary waves 

For p slightly larger than p t , the homogeneous solutions w = and w = wi are 
unstable, and one should look for the onset of spatio-temporal structures rather than 
purely stationary states. In this respect, one may be guided by the observations made 
in numerical simulations [141 115] . where 'stripes' of higher density moving over a low 
density background have been reported. Such structures are rather similar to soliton 
solutions that have been observed in many different physical contexts [2Tj . 

5.3.1. Stationary hydrodynamic equation in a moving frame. Let us now look for 
possible soliton solutions of the hydrodynamic equations (jlip and (|27|) . To this aim, 
we assume for p and w the following "propagative" form, with propagation velocity 
c > 0, along an arbitrary axis x of unit vector e: 

p(r,t) = R(x - ct), w(r,t) =W(x- ct)e, (73) 

with £ = x — ct and W(() > 0. Using Eqs. (fTT) , one finds the simple relation 
R' = W'/c, leading to: 

R(0 = lw(0+P*, (74) 

where p* is up to now an arbitrary constant density. In the following, we consider 
velocity profiles that vanish for £ — > ±oo, so that p* = lim^ioo R{(). Inserting this 
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form in Eq. Q27[) . one can eliminate i?(C) and obtain the following ordinary differential 
equation for W(£), also taking into account the density dependence of the transport 
coefficients [§]: 

W" = -(ao - a x W - a 2 W')W - b{W - b 2 W 2 - b 3 W 3 . (75) 
The different coefficients in Eq. (|75| read 

ao = (c- V f\ (Di + D 2 p*) (76) 



ai=7 + ^2^-lj (77) 

a 2 = — r (78) 

c(D 1 + D 2 p*) V ' 

h = p J '(p*-p t )(D 1 +D 2 p*) (79) 

b 2 = ^[D 1 +D 2 (2p* -p t )] (80) 
c 

h = — £ (81) 



with 



D 1 = ^(l-e-^) (82) 

D2 =3^U +e )' < 83 > 

and 7 = 7/V, £ = £,/v. As often in the study of solitons [21], Eq. (|75p may be 
reinterpreted as the equation of motion of a fictive particle with position W at time 
Q. Here, this virtual particle has a unit mass, and moves in a potential 

= h -±W 2 + ^W 3 + h -^W\ (84) 

with a non-linear friction force — (ao — a\W — a 2 W')W . Depending on the sign of the 
effective friction coefficient (a — a{W — a 2 W), the friction force may either dissipate 
or supply energy to the particle. Note that this friction term breaks the symmetry 
£ — > — C, so that the resulting momentum profile cannot be symmetric. 



5.3.2. Numerical integration of the velocity and density profiles. To find a solution 
for W(£), we integrate numerically Eq. (|T5[) for given values of the parameters ai and 
bi. The following constraints are imposed to the solution: W(Q should be positive for 
all values of £, and W(C) should go to for £ — > ±00. Hence for large values of 
W(Q should be small, and should satisfy, to a good accuracy, the linearized version 
of Eq. (ITS]) , namely: 

W" + aoW + hW = 0. (85) 

This equation has two exponential solutions W±((^) = A± exp(r±£), with: 

r± = ~ (-ao ± sJal-Ab^j . (86) 

§ We however neglect the density dependence of the ratio v 1 /u, as it would lead to terms of higher 
order than that retained in our hydrodynamic description. 
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Figure 7. (a) Velocity profile v(x, t) = V(£), with £ = x — ct, for p* = 0.06 and 
a = (to = 0.31 (dashed line), 0.32 (dot-dashed line) and 0.33 (full line). Inset: 
propagation velocity c as a function of <r. (b) Density profile p(x, t) = R(C) for the 
same values of the parameters. Horizontal dotted lines correspond to the density 
pt for cr = ctq = 0.31, 0.32 and 0.33 (bottom to top). Other parameters: A = 0.5, 
do = 0.5 and vq = 1. 



For W(C) to be positive, one needs that the roots r± be real, which implies Oq — 4b\ > 0. 
Further, for W(£) to vanish both for £ — * — oo and ( — > +oo, one should have both an 
increasing and a decreasing solution for Eq. (|85p . namely r + > and r_ < 0, which 
corresponds to fri < 0. 

The free parameters in Eq. (|75j) are c and p* (this point will be briefly discussed in 
Section 15.41 in connection with numerical results). The noises a and cto are external 
control parameters. The overall density p is computed afterwards from the profile 
R(C)- Assuming that we are in the low noise region of parameter space a < cr£°, then 
[exp(— ct 2 /2) — 2/3] > and the condition b\ < implies p* < pt- In addition, as 
the trajectory of the particle starts and ends at the same position W = with zero 
velocity (W — 0), its energy is the same, which means that the friction force has 
to dissipate energy on some part of the trajectory and to supply energy otherwise. 
Assuming oq > implies c — Vq/(2c) > 0, that is c > vq/\/2. On the other hand, one 
intuitively expects c to be smaller than the microscopic velocity vq of the particles. 

The numerical procedure we implement is the following. Choosing a given value 
for p* and for c, we start at £ = ( < (|Co| > 1), with a small value W(£o) = Wo « 1. 
and with a derivative M^'(Co) = r+Wo- This choice of initial conditions ensures that 
we select a solution with an exponential tail W(() = A + exp(r+£) for £ < (b- Then we 
integrate numerically the differential equation for £ > Co, until reaching large enough 
positive values of £. At this stage, two behaviours may appear. One should first 
notice that for b% < (and at least if &3 does not take a large negative value) the 
effective potential has a local maximum in W = and a local minimum at a 

value W = Wmin- Then, if the dissipated energy is larger than the injected energy, the 
particle ends up at W m i n for ( — > oo, yielding a profile W(Q that does not fulfil the 
condition required. In the opposite case, if energy injection dominates, the particle 
crosses the local maximum at W — and goes to negative values. It is only in the 
marginal case where dissipation exactly compensates injection that the correct profile 
W(C) is found. As friction is mainly controled by the parameter c, we keep p* fixed and 
perform a loop over the value of c in order to converge toward the marginal solution. 
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Figure 8. Solitons in the numerical model, (a) Instantaneous snapshot, the band 
is moving south-west; lengths are scaled by do. (b) Example of trajectory in the 
direction of the averaged velocity, (c) Mean profiles along the direction of the main 
motion. We plot the reduced dimensionless density R r = ((p(x — ct)) — p sat )dgc/«o 
(dotted line) and the dimensionless momentum Wi = (w(x—ct))d'^ ) /vQ (plain line), 
both being time-averaged in the comoving frame of the soliton. (d) and (e) Same 
data as (c) on semi-log scales, emphasizing the exponential decay. The scales 
are identical on vertical axes, but different on abscissas. Parameter values are 
p = 2 , <r = 0.163, L = 4096; the other ones correspond to set VII in Tablc[l] 



Note however that if 63 < 0, <&(W) — > —00 when W — > +00, so that one should also 
take care that the particle does not "escape" to large positive values of W. 

Using the above procedure, we obtain a family of profiles W(() with three control 
parameters, namely the "background" density p* and the noises a and <jq. In the 
following, we restrict ourselves to the case 00 = a . The density profile is computed 
from the relation R(Q = P* + W(()/c, and the velocity profile V(Q is obtained from 
the momentum profile W(Q through V(Q = W(£)/R(£). Examples of such velocity 
profiles are presented in Fig. [3 for different values of er and for a given value of p* . 

A remaining open issue is the stability of these solitary waves with respect to small 
perturbations. A formal stability analysis like the one performed for the homogeneous 
state of motion is a difficult task here, and we leave this question for future work. 

5.4- Solitary waves in the agent-based model 

We now compare the solitary waves computed in the hydrodynamic equations with the 
travelling stripes observed in direct numerical simulations of the agent-based model 
(see Fig. (Hta))- We focus again on the case 00 = er. The stripped structures are 
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Figure 9. Solitons in numerical model, (a) Density of saturated vapour for 
different p (p = 2~ 2 , o, 5 _1 , ■. 2~ 3 , A, 10~\ T and 2 -4 , *). Inset: finite size 
effects on soliton, p = 2~ 4 , L = f024 X and L = 2048 +. The dashed line 
marks the value of the global density (p = 2 — 5 ). The dotted lines underline the 
threshold of the collective motion at and of the homogeneous moving population 
er r .(d) Speed of the solitons (same parameters). The other parameters are the 
ones of set VII (see table [TJ. 



composed by several localised, randomly spaced bands. They are not part of a 
regular pattern, nor a wave train |15j . They are all moving along the direction of 
the main motion, although during the transient period they can pass through each 
other with only few interactions. The space between two bands is filled with particles 
moving independently (the hydrodynamic momentum vanishes), and homogeneously 
(the density is constant). In analogy to the liquid-gas coexistence, we denote this state 
as the saturating vapour. 

We observe that the bands move at a constant speed, at least on the duration 
necessary to travel through the system size (Fig. |S[b)). From the trajectories, we 
measured the velocity c of the solitons. On the density profiles, we extracted the 
value p sat of the density outside the peak. If these structures are only propagative and 
if the continuity equation is valid at a coarse-grained level in the agent based model 
(which is expected from mass conservation) , the density and momentum profiles should 
be related by W = c(R — p sat ), as in Section IQ1 Plotting on Figure [S^c) both the 
reduced density c(R — p sat ) and the momentum W, we observe that both curves match 
perfectly, confirming the propagative nature of this stripped pattern. 

These solitary waves are quite similar to the soliton we found in the hydrodynamic 
equations (see section f5.3|) . with in particular an exponential decay of the momentum 
profile on both sides (Fig. [SJd-e)). However, the asymmetry of the profile is much 
more pronounced than in the analytical model: the exponential decay is much steeper 
in front of the profile than in the rear part. 

We now study how the two main characteristics of the solitary waves, namely the 
velocity c and the density p sat , vary with the control parameters p and a. Since we 
perform a numerical study, we need to be aware of finite size effects. Plotting the 
density of saturating vapour versus noise for a given density but for different sizes, we 
can make three observations (inset of Fig. [9]). First, the system size hardly changes 
p sat provided that the noise amplitude remains near the threshold. Moreover, at a 
lower noise amplitude, the density p sat increases and become sensitive to the system 
size. Lastly, there is a noise c r below which we cannot observe solitons anymore (see 
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also [15]) and the system becomes homogeneous at a coarse-grained level. The result 
is qualitatively consistent with the restabilization of the homogeneous flow described 
in Sections 15. II and 15.21 

The study of the very low noise amplitude region of the phase diagram is an 
ongoing work. So we mainly focused in the present article on the region relatively 
close to the transition to collective motion. For different global densities, both p sa± 
and c fall onto the same curve when plotted as a function of a, as shown on Fig. [9]^ a) 
and (b). Therefore, once the noise amplitude a is given, the characteristics (c,p sat ) 
of the solitary waves are determined, and the number of solitary waves is adjusted by 
the dynamics in order to match the global density of the system. 

This is a major difference with the solitary waves obtained from the hydrodynamic 
equations in section 15.31 These solitary waves depend on two control parameters, 
namely the noise amplitude a and the density at infinity p* . Hence there is a priori 
no way to determine the number of solitons in a large but finite system with a given 
density. At a heuristic level, we might guess that the solitary waves may be stable 
only for some specific values of p*, which would give a selection mechanism for the 
density p sat . Such a mechanism would make the connection between the analytical 
and numerical models clearer, but we presently have no clue to confirm this tentative 
scenario. Obviously, further studies of the dynamics of the solitary waves in the 
context of the hydrodynamic equations are needed. 

6. Conclusion 

In summary, we have derived in this article hydrodynamic equations for a model of 
self-propelled particles with binary interactions, in the regime of low hydrodynamic 
velocity. We also compared the results of the hydrodynamic description to the 
numerical simulations of a standard agent-based model. In the analytical model, 
the homogeneous state with zero velocity is a stationary solution for any values of the 
microscopic parameters (the noise amplitude and the overall density), but this state 
is linearly unstable for a reduced density p greater than a transition density pt (a) , or 
equivalently, for a noise smaller than a transition value <7t(p). 

When the zero velocity solution is unstable, another homogeneous state, with a 
nonzero hydrodynamic velocity, appears. This state is linearly stable with respect to 
spatially homogeneous perturbations. However, close to the transition line <7 t (p), this 
state turns out to be linearly unstable with respect to finite wavelength perturbations. 
As the validity of the hydrodynamic equations is, strictly speaking, restricted to the 
vicinity of the transition line, we also studied the stability of the homogeneous state 
of motion directly from the Boltzmann equation. We found that, far enough from 
the transition line, the homogeneous motion becomes linearly stable. Interestingly, 
this restabilization phenomenon is also qualitatively observed in the hydrodynamic 
equations, although this regime is beyond their domain of validity. All these results 
agree semi-quantitatively with the numerical simulations of the agent based model. 

When the homogeneous state of motion is unstable, more complex spatio- 
temporal structures should appear. A candidate for such structure is the solitary 
waves we obtained from the hydrodynamic equations. These solitary waves resemble 
the moving stripes observed in the numerical agent-based model, apart from the 
asymmetry which is more pronounced in the latter. A word of caution is however 
needed here, as on the one hand the solitary waves have a finite amplitude, so that 
the hydrodynamic equations might not be valid, and more importantly, their stability 
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has not been tested yet. On the basis of the numerical simulations of the agent-based 
model, it is however likely that these solitary waves should be stable at least in a given 
region of the phase diagram. 

As for future work, it would be interesting to investigate the stability of the 
solitary waves, and to look for possible "multi-soliton" solutions, in case the stability 
would be confirmed. Specifically, it would be interesting to be able to determine the 
number of solitons, their celerity and the background density as a function of the global 
density (for a finite volume) and of the noise amplitude, if such a relation exists, as 
suggested by the numerical simulations of the agent-based model. 
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Appendix A. Agent-based model 

Appendix A.l. Looking at the model further 

The numerical system we looked at is very similar to the one defined by Vicsek 
et al [12J. This is a very minimal model, easy to implement. In contrast, a real 
direct simulation would have been coded following a molecular dynamics algorithm, 
which would have cost much more cpu time than our Monte Carlo-like program. 
The numerical choice is also related to the fact that collective motion of self- 
propelled particles has been mainly studied in this framework during the last ten 
years [21[21[ll[I3[53[Il[151[23[Il[2Hl[2Sl[3ni- Thus we would like to take profit 
from this large background and the knowledge of the system we already got. 

To fully understand the results presented in this paper, we must explain the 
differences between the numerical system we used and a direct simulation. In what 
we have done, collisions are computed at fixed time step. So every other collision that 
could have occurred within At is neglected. On the other hand, collisions can involved 
many individuals. Another implication of the discrete time step is that decreasing the 
time step increases the collision frequency. Then the noise does not act on the system 
with the same manner for two different time steps. Hence, in its present formulation, 
the agent-based model is not a discretized version of a continuous time model. To 
reach this goal, the noise amplitude should be renormalized in some way with the 
time step. 

The balance of the above different effects is difficult to imagine a priori. We 
do not expect any quantitative matching between the theory we developed and the 
simulations we presented. But we still want to test the robustness of the predictions 
made for large system sizes. 

We must also emphasize that some studies in the literature were aimed at giving 
an exact continuous theory of Vicsek's model [3lJ [32] [33]. Up to now, this difficult 
problem has been dealt in the framework of perturbative theories at a first order in 
speed differences. In addition, the role of the noise is not properly taken into account 
in these studies: it is either ignored [31] , or described by a phenomenological diffusive 
term 32 . Finally, the transport coefficients of the hydrodynamic equation do not 
contain any dependence on the microscopic parameters of the model. 
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Figure Al. Phase transition in numerical simulation and finite size effects. We 
plotted (a) the Binder cumulant K and (b) the averaged order parameter (ip) vs 
noise rms-value, for three different sizes. In the inset, we show the histogram of 
the order parameter ip at the transition point for a system size L = 1024. On 
figure (a), we emphasized the depth of the well Ait: approximation of the errors 
in determining the transition point. Parameters are the ones of configuration n° 
III with p = 1/16 or p = 1/2. 



Appendix A. 2. Phase transition 



As in usual versions of self-propelled particles systems, the behaviour of the system 
roughly falls into two different categories. Either there is no collective motion: every 
particle move randomly without clear correlation with its neighbours; or there is a 
non-zero global velocity in an arbitrary direction. 

Since the analogy with magnetic systems is quite obvious, the habits is to consider 
the equivalent averaged magnetization of our system, namely the global normalized 
velocity ip*: 



1 



Nv, 



N 

£vj 



(A.l) 



considered as an order parameter. To determine the characteristics of the phase 
transition, we study the statistical properties of the order parameter (/?', considering 
its mean (ip), its variance x ancl its Binder cumulant K |34j : 



{<P) 



T 



t=i 



K = 1- 



L 2 ((^) - 



(A.2) 

(^> 2 ) , (A.3) 

W (A ' 4) 
The brackets (. . .) indicate an averaging over time. The duration of the simulation 
has to be large to inhibit memory effects. Ideally, the correlation time for each 
set of parameters (p, a, o~o, vq, do, A) should be computed from the auto-correlation 
function [35]. However, this is a tantamount tasl|jj|. Practically, in order to have a 
rough approximation of the correlation time, we measured the transition time from 



The cumulative consumed cpu time already reaches fifty years. 
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the initial condition to the stationary state. Then we performed averaging on time 
which are hundred times greater than that transition time. 

For all sets of parameters I to VII (Table [T]), we observed that the system exhibits 
a phase transition from a non-moving to a globally moving population when decreasing 
the noise amplitude at a fixed density. At small enough size L, all statistical variables 
((ip) , x, K) remain continuous, while a singular point appears when the system is larger 
than a typical size L t (see Fig. lAlf a) and (b), as well as Refs. [HI HI]). 

The main observations are the following: the order parameter curve exhibits a 
jump (Fig. lAlf b)). the variance is delta-peaked (not shown here), the Binder cumulant 
has a minimum (Fig. lAlf a)) which goes to larger negative values when the system 
size is increased, and the histogram of the order parameter is bimodal (see inset of 
figure lAlT b)). All these sign plead in favour of a first-order phase transition. 

It is now well known that a finite size system exhibits a rounded transition, at 
equilibrium [36, 37 or far from the equilibrium 38, 391. The scaling laws are normally 
sufficient to detect the order of the transition. In our case, the finite size scaling 
laws correspond to a continuous transition below L t (40] [41] , and to a discontinuous 
transition above L t [15j . 

To estimate the transition point, we measured the location where the Binder 
cumulant minimum becomes negative. We neglected the finite size effects at higher 
size. We determined the error bars on that location as the width of the well (see 
Fig.EIa)). 



Appendix B. Stability against arbitrary perturbations 



In this appendix, we study within the framework of the hydrodynamic equations the 
stability of the stationary homogeneous flow. Starting from Eqs. (j41]) and (|42|) . we 
consider the case where wo = wi ^ 0, solution of Eq. (|33[). The rotational symmetry 
is broken when the collective motion appears, and we take e|| = wi/|wi| as a first 
vector of the geometrical basis. Then we define two angles $i and $2 between e|| and 
the directions of 5wq and q respectively. We denote as ej_ the unit vector orthogonal 
to e||, and such that (e||,ej_) form a direct basis. 

From Eqs. ([44]) and (|45|) . we project the resulting vectorial equation onto ej| and 
ex, and we eliminate the ratio Swo/Spo from the continuity equation, yielding: 



s is - 



i"/w±q cos $2 



1 



|^gcos(^i 



vq \ cost^i = 
0a) 



(B.l) 



isKWi cos??] 



q cos $2 



[2s£,wl cos $1+ i (^' — £'wl + sk) qw\ cos($i — $2)] , 



[s + i"/Wiqcos$2 + vq 2 



--w gcos(i?i 



sin$i = 
0a) 



(B.2) 



isKWi cost?] 



q sin -#2 



where q and wi are real and positive. These are two polynomial equations that we 
will study at a given point (p, a) of the phase diagram, for a set of physical variables 
(do,A,t>o) an d f° r different pairs (^1,^2)- For all fixed parameters, the solutions of 
those equations will be a discrete number of sets (5, s). 

First, one can check that this set of equations is invariant when Sw is rotated 
with an angle of it (t?i — > i?i +7r). Note also that every real term depends on an even 
power of q. So one can expect that the real part of the growth rate 3? [s] only depends 
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Figure Bl. Stability against inhomogeneous perturbations, solutions of 
Eqs pO)) and lHO) so that S[<j] = 0. (a) (#1,02), with a = ct = 0.5, p = 0.22 
(dot-dashed line), 0.30 (dashed), 0.40 (full line), (b) 5R[s] vs q, same parameters 
as (a), Sf[s] increases when p increases. 



on even powers of q, and that Sft[s] remains invariant when $2 is changed into $2 + tt. 
That is why we will study Eqs. fBTT]) and (|R2|) for (0i,i? 2 ) € [0, tt[x [0, tt[. 

A third property arises clearly when we introduce the expressions ([25)1 - (l3"2j) ) of the 
transport coefficients in Eqs. (|B.1[) and (|B.2[) : the wavenumber q appears only through 
the product qBado, meaning that the solutions q are proportional to l/Bdo. As 
already mentioned, the framework of the kinetic approach implies that B is large, and 
therefore it implies that we are studying long wavelength perturbations. This analysis 
also shows that the growth rate depends only on the dimensionless control parameter 
p, the noise amplitudes a and ctq, and the self-diffusion rate A which gives the proper 
unit to s. Let us also mention that a trivial solution of the system of equations is 
s = for q = 0. This solution is actually an artefact of the calculation procedure 
(namely a multiplication by s), as it is not a solution of the original equations ()41|) 
and P^]) . Hence we will not consider this extra solution in the following. 

For some parameters , 1^2)5 one or several terms can vanish, and the degree of 
the polynomials may decrease. We will first study the general equations and those 
particular cases will be considered later. If we combine ( (IB. lft x sin rtB.2|) x cos 
we get a linear equation in s : 



qvg sin(i?i - # 2 ) + 2iw 1 (// - £'wf) siwdi 



2£wi cos #1 sin d\ + iqn sin d 2 



Indeed, we verify that 3? [s] is an even function of q. Now we can replace s by its 
expression in Eq. (jB.2|) . The resulting equation is a third degree polynomial, that can 
be formally written as: 

d 3 q 3 + id 2 q 2 + d x q + id Q = 0, (B.4) 

where the coefficients di are real functions of (p,a,ao, B) and of , i?2 ) • The last 
three coefficients are rather difficult to manipulate. But this equation can be easily 
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solved, using Cardano's method for instance. In the case where 
d 3 = sin$i sin $ 2 sin(#i — $ 2 ) cos($i — $2) 



(B.5) 



does not vanish, we compute the solutions. The resulting values q are complex 
numbers, so that they cannot correspond to physical solutions. Yet, for some sets 
of angles ($1,^2), the solutions for q are real. We are interested only in these modes. 
Determining the angles (i9i,i?a) f° r which q is real, we then compute the growth 
rate 3i[s+] using Eq. (|B.3[) . There are four different branches (see Figs. IBlf a)). whose 
lengths increase when the control parameter is chosen deeper in the collectively moving 
phase (i.e. at low a or at high p). For all sets of parameters for which we have 
computed the growth rate, its real part remains negative (Fig. IBlf b)). Thus the 
homogeneous moving phase is stable against finite wavelength perturbations in the 
general case. 

The above calculation relies on the assumption that d% 7^ 0. This assumption is 
not valid in either of the four following cases: 

• a longitudinal instability (sin$i = 0), 

• a wave vector q colinear to the direction of the main motion (sin$ 2 = 0), 

• a perturbation (5wo colinear to the wave vector q (sin(i?i — $2) = 0). 

• a perturbation (5wo perpendicular to the wavevector q (cos($i — $2) =0). 

We first consider the study of stability under a longitudinal perturbation: Wj and 
5w are colinear. Then the Eq. (|B.2|) vanishes in two cases: 



Replacing s by the first expression in equation (|B.1|) . we can show that there is no 
authorised mode, in other words ^s[q] 7^ 0. So, Eq. (|B.2[) vanishes only for #1 = $2 = 0. 
The corresponding stability analysis is presented in details in Section 14.2.31 

For any of the last three cases, we solve equations IB. 21 and IB. 21 and we find that 
either there is no authorized mode (q is complex), or 9?[s+] < 0. Thus none of those 
cases is related to an unstable mode. 

To sum up, this study of the stability of the homogeneous stationary moving 
phase shows that the longitudinal direction is the only mode which can be unstable. 
This result is consistent with the observations made in numerical simulations [14[ I15j . 
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